log using "$rootpath/log/Fig4.log", replace
//Figure 4 final run 

use "$rootpath/processed_data/Estimation_Data.dta", clear
drop if key_estimation_sample!=1
sort facility_id year
//Note the quartiles in the other spec are 0.475 and 0.24
egen num_obs=nvals(year), by(facility_id)
sort year
gen temp_damage=log_damage55 if num_obs==10
egen mean_damage=mean(temp_damage), by(year)
gen temp_damagelow=log_damage55 if num_obs==10&avg_old_code_prob<0.23
egen mean_damagelow=mean(temp_damagelow), by(year)
label var mean_damagelow "Plants with Low Treatment" 
gen temp_damagehigh=log_damage55 if num_obs==10&avg_old_code_prob>0.53
egen mean_damagehigh=mean(temp_damagehigh), by(year)
label var mean_damagehigh "Plants with High Treatment"
label var mean_damage "All Plants"
label var year "Year"
//summarize mean_damage  if year==2014
//summarize mean_damage  if year>2014
// display ((1.28042 - 1.325634)/(1.325634))
//-.03410745
line mean_damagelow year if num_obs==10&avg_old_code_prob<0.23,  lpattern(longdash) ytitle("Inverse Hyperbolic Sine (damages 100,000 USD)",height(7)) ||  line mean_damagehigh year if num_obs==10&avg_old_code_prob>0.53, lpattern(shortdash)|| line mean_damage year if num_obs==10, yscale(r(1 2))  xlabel(2010 "2010" 2011 "2011" 2012 "2012" 2013 "2013" 2014 "2014" 2015 "2015" 2016 "2016" 2017 "2017" 2018 "2018" 2019 "2019")  title("Plant Damages") subtitle(" ") xline(2014) legend(bmargin(text)) graphregion(color(white))
graph export "$rootpath/table_graph/Figure_FourB.pdf", replace


use "$rootpath/processed_data/Estimation_Data.dta", clear
drop if key_estimation_sample!=1
sort facility_id year
//Note the quartiles in the other spec are 0.475 and 0.24
egen num_obs=nvals(year), by(facility_id)
sort year
replace tot_damage=tot_damage/1000000
gen temp_damage=tot_damage if num_obs==10
egen mean_damage=mean(temp_damage), by(year)
gen temp_damagelow=tot_damage if num_obs==10&avg_old_code_prob<0.23
egen mean_damagelow=mean(temp_damagelow), by(year)
label var mean_damagelow "Plants with Low Treatment" 
gen temp_damagehigh=tot_damage if num_obs==10&avg_old_code_prob>0.53
egen mean_damagehigh=mean(temp_damagehigh), by(year)
label var mean_damagehigh "Plants with High Treatment"
label var mean_damage "All Plants"
label var year "Year"
line mean_damagelow year if num_obs==10&avg_old_code_prob<0.23,  lpattern(longdash) ytitle("Damages (millions USD)",height(7)) ||  line mean_damagehigh year if num_obs==10&avg_old_code_prob>0.53, lpattern(shortdash)|| line mean_damage year if num_obs==10, yscale(r(2 6))  xlabel(2010 "2010" 2011 "2011" 2012 "2012" 2013 "2013" 2014 "2014" 2015 "2015" 2016 "2016" 2017 "2017" 2018 "2018" 2019 "2019")  title("Plant Damages") subtitle(" ") xline(2014) legend(bmargin(text)) graphregion(color(white))
graph export "$rootpath/table_graph/Figure_FourA.pdf", replace

use "$rootpath/processed_data/Estimation_Data.dta", clear
drop if key_estimation_sample!=1
sort facility_id year
//Note the quartiles in the other spec are 0.475 and 0.24
egen num_obs=nvals(year), by(facility_id)
sort year
//replace tot_damage=tot_damage/1000000
gen temp_damage=tot_emissions if num_obs==10
egen mean_damage=mean(temp_damage), by(year)
gen temp_damagelow=tot_emissions if num_obs==10&avg_old_code_prob<0.23
egen mean_damagelow=mean(temp_damagelow), by(year)
label var mean_damagelow "Plants with Low Treatment" 
gen temp_damagehigh=tot_emissions if num_obs==10&avg_old_code_prob>0.53
egen mean_damagehigh=mean(temp_damagehigh), by(year)
label var mean_damagehigh "Plants with High Treatment"
label var mean_damage "All Plants"
label var year "Year"
line mean_damagelow year if num_obs==10&avg_old_code_prob<0.23,  lpattern(longdash) ytitle("Criteria Emissions (tons)",height(7)) ||  line mean_damagehigh year if num_obs==10&avg_old_code_prob>0.53, lpattern(shortdash)|| line mean_damage year if num_obs==10, yscale(r(0 200))  xlabel(2010 "2010" 2011 "2011" 2012 "2012" 2013 "2013" 2014 "2014" 2015 "2015" 2016 "2016" 2017 "2017" 2018 "2018" 2019 "2019")  title("Plant Emissions") subtitle(" ") xline(2014) legend(bmargin(text)) graphregion(color(white))
graph export "$rootpath/table_graph/Figure_FourC.pdf", replace


use "$rootpath/processed_data/Estimation_Data.dta", clear
drop if key_estimation_sample!=1
sort facility_id year
//Note the quartiles in the other spec are 0.475 and 0.24
egen num_obs=nvals(year), by(facility_id)
sort year
gen temp_damage=log_emissions4 if num_obs==10
egen mean_damage=mean(temp_damage), by(year)
gen temp_damagelow=log_emissions4 if num_obs==10&avg_old_code_prob<0.23
egen mean_damagelow=mean(temp_damagelow), by(year)
label var mean_damagelow "Plants with Low Treatment" 
gen temp_damagehigh=log_emissions4 if num_obs==10&avg_old_code_prob>0.53
egen mean_damagehigh=mean(temp_damagehigh), by(year)
label var mean_damagehigh "Plants with High Treatment"
label var mean_damage "All Plants"
label var year "Year"
line mean_damagelow year if num_obs==10&avg_old_code_prob<0.23,  lpattern(longdash) ytitle("Inverse Hyperbolic Sine (emission tons)",height(7)) ||  line mean_damagehigh year if num_obs==10&avg_old_code_prob>0.53, lpattern(shortdash)||  line mean_damage year if num_obs==10, yscale(r(1 2))  xlabel(2010 "2010" 2011 "2011" 2012 "2012" 2013 "2013" 2014 "2014" 2015 "2015" 2016 "2016" 2017 "2017" 2018 "2018" 2019 "2019")  title("Plant Emissions") subtitle(" ") xline(2014) legend(bmargin(text)) graphregion(color(white))
graph export "$rootpath/table_graph/Figure_FourD.pdf", replace

log close
